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Circadian oscillation provides selection advantages through S3mchronization 
to the daylight cycle. However, a reliable clock must be designed through 
two conflicting properties: entrainability to synchronize internal time with 
periodic stimuli such as sunlight, and regularity to oscillate with a precise 
period. These two aspects do not easily coexist, because better entrainability 
favours higher sensitivity which may sacrifice regularity. To investigate condi- 
tions for satisfying the two properties, we analytically calculated the optimal 
phase -response curve with a variational method. Our results indicate an 
existence of a dead zone, i.e. a time period during which input stimuli neither 
advance nor delay the clock. A dead zone appears only when input stimuli 
obey the time course of actual solar radiation, but a simple sine curve 
cannot yield a dead zone. Our calculation demonstrates that every circadian 
clock with a dead zone is optimally adapted to the daylight cycle. 
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1. Introduction 

Circadian oscillators are prevalent in organisms from bacteria to humans and 
serve to synchronize bodies with the environmental 24 h cycle [1,2]. Although 
the molecular implementation of oscillation is species-specific [3-6], every circa- 
dian clocks satisfies two requirements to achieve reliable synchronization to the 
environment: entrainability to synchronize internal time with periodic stimuli 
and regularity to oscillate with a precise period. Circadian clocks are acquired 
through evolution independently in bacteria, fungi, plants and animals [7]. None- 
theless, entrainability and regularity constitute major characteristics conserved in 
all circadian clocks [6], which strongly suggest that these two properties are essen- 
tial for survival. A main source of interference with regularity is discreteness of 
molecular species, i.e. molecular noise [8-13]. Many studies have analysed the 
resistance mechanisms of circadian oscillators against noise [14-17]. Regarding 
entrainability, circadian clocks synchronize their internal time with the environ- 
mental cycle via sunlight, and its effect depends on the wavelength or fluence, 
as well as on the phase of the stimulation. However, entrainability and regularity 
are conflicting factors, because circadian clocks with better entrainability are sen- 
sitive not only to the periodic light stimuli, but also to the molecular noise which 
interferes with regularity. 

Because both regularity and entrainability are important adaptive values, we 
expect actual circadian oscillators to optimally satisfy these two factors (figure 1). 
Here, we investigate the optimal phase -response curve (PRC), which is both 
entrainable and regular, in the phase oscillator model [18] by using the Euler- 
Lagrange variational method. Our main finding is the inherent existence of a 
dead zone in the PRC: optimality is achieved only when the PRCs have a time 
period during which light stimuli neither advance nor delay the clock (figure la). 
In other words, a PRC with a dead zone (figure Tm) is better adapted than those 
without a dead zone (figure 2h). This result is intriguing, because a dead zone, 
with which oscillators tend to be unaffected by stimuli (i.e. lower entrainability), 
achieves better entrainability. We also tested this with two types of input stimuli: 
a solar radiation-type input that simulated the time course of solar radiation 
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Figure 1. Illustrative relation between two trade-off properties: entrainability 
and regularity. There is an infeasible region with respect to entrainability and 
regularity (coloured area), inside which no clocks can be implemented. Actual 
circadian clocks are considered to optimally satisfy them and such optimal 
clocks lie on the edge between feasible and infeasible regions (thick- 
dashed line). 



intensity (cf. equation (2.24) and figure 4fl) and a simple 
sinusoidal input (sine curve). Surprisingly, the dead zone in 
the optimal PRC emerges only for the solar radiation-t3rpe 
input, not for the sinusoidal input. Many experimental studies 
reported the existence of a dead zone in various species (figure 
2c 4 show experimentally observed PRCs of fruitfly [19] and 
mouse [20], respectively). Our results indicate that circadian 
oscillators in various species have adapted to solar radiation 
for reliable synchronization. 



2. Models and methods 
2.1. Phase oscillator model 

Circadian oscillators basically comprise interaction between 
mRNAs and proteins, whose dynamics can be modelled by 
differential equations. A circadian oscillator of N molecular 
species can be represented by 

dxi 
dF 



:F,{x) (z = l,2,...,N), 



(2.1) 



where the N-dimensional vector x = (xi, X2, . . . ,X]^) denotes 
the concentration of molecular species (mRNAs or proteins). 
The effect of noise on genetic oscillators has been a subject of 
considerable interest, and noise-resistant mechanisms have 
been extensively studied [14-17,21-23]. In general, the 
dynamics of the zth molecular concentration in a circadian 
oscillator subjected to molecular noise is described by the 
following Langevin equation (Stratonovich interpretation): 

dxi 



dt 



= F,{x-p) + Q,{x)m. 



(2.2) 



where Qi{x) is an arbitrary function representing the multi- 
plicative terms of the noise, ^i{t) is white Gaussian noise 
with the correlation {^i{t)^j{t')) = 28ij8{t - t') (a bracket ( • ) 
denotes expectation), and p is a model parameter. 

Circadian oscillators synchronize to environmental cycles 
by responding to a periodic input signal (light stimuli). We let 
p in equation (2.2) be stimulated by the input signal: for 
example, p can be the degradation rate (for the sake of simpli- 
city, we consider that the input signal affects only one 
parameter). We use equation (2.2) for calculating regularity 
and entrainability of circadian oscillators. 



2.2. Definition of regularity 

Because the circadian oscillator of equation (2.2) is subjected 
to noise, its period varies cycle to cycle. We use the term 
regularity for the period variance of the oscillation (higher 
regularity corresponds to smaller period variance). Let us 
first consider the case without input signals (i.e. p is cons- 
tant). As equation (2.1) exhibits periodic oscillation, we can 
naturally define the phase (/> G [0,2 tt) on equation (2.1) by 



dcj) 

'dt' 



a 



(2.3) 



where il = 27r/T is the angular frequency of the oscillation 
(T is a period of the oscillation). The phase (/> in equation (2.3) 
is defined only on a closed orbit of the unperturbed limit- 
cycle oscillation. However, we can expand the definition into 
the entire x = {xi, X2, . . . ,x^) space, where the equiphase surface 
is referred to as the isochron X( (/>) (figure 3a). By using standard 
stochastic phase reduction [18], equation (2.2) can be trans- 
formed into the following Langevin equation with respect to 
the phase variable cj) (Stratonovich interpretation): 



d^ 

'dt' 



/2 + £^-(c/>)QK0)ft(O, 



(2.4) 



where IJ((/)) = (Ui( (/>),... ,U7^((/>)) is an infinitesimal PRC 
(iPRC) U((/)) = ^x4>\x=xi^c{^)' abbreviated Q/(xlc(</>)) as 

Qi{4>). iPRC Ui{4>) quantifies the extent of phase advance or 
delay when perturbed along an Xi coordinate direction at 
phase (/). The N-dimensional vector Xi^c{4*) denotes a point on 
the limit-cycle trajectory at phase (/>, where LC stands for limit 
cycle. The value of iPRC Ui{4>) is calculated as a solution of 
an adjoint equation [24] or as the set of eigenvectors of a mono- 
dromy matrix in the Floquet theory [18] for arbitrary oscillators. 
Let P{4>;t) be the probability density function of 4> at time t. 
From equation (2.4), the Fokker-Planck equation (FPE) [25] of 
P{4>;t) is given by 



aP(c/>;0 



dt 



-^(/2 + ^(c/>))+^e(c/>) ^P(c/>;0, 



(2.5) 



where 



and 



'd4, 



Ui{mi{<i>) (2.6) 



(2.7) 



Introducing a slow variable ip= 4>- fit, the FPE of the 
probability density function n{(p; t) = P{4>= + fit; t) is 
given by 



dt 



n(cp;t) 



- ^^((P + nt) + + m }n{cp; t). 



(2.8) 



With sufficiently weak noise, n{cp;t) is a slowly fluctuat- 
ing function of t. In such cases, T{cp + fit) and G{cp + fit) 
fluctuate much faster than n{(p;t), thus these two terms can 
be averaged for one period while keeping n{cp;t) constant 
(phase averaging). In other words, we separate time scales 
between T((p -\- fit), Q{(p -\- fit) and n{(p;t). By phase aver- 
aging, T{ip + fit) vanishes because of the periodicity (use 




Figure 2. Illustrations of typical PRCs (a) with and [b) without a dead zone. Experimentally observed PRCs as a function of time in (c) fruitfly [Drosophila) [19] and 
id) mouse (Mus) [20] with light pulses (circles) and their trigonometric fitting curves (solid line). Shaded and non-shaded regions indicate subjective night and day, 
respectively. 
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Figure 3. (a) Illustration of the isochron Xicj)), the solid and dashed lines describing a limit-cycle trajectory and its isochron drawn at intervals of 7t/6, respectively, (b) 
Relation between the phase variance V^/, and the period variance Vt in Langevin equation (2.2) (the solid lines represent trajectories of the Langevin equation). V^^ is the 
variance of the phase cj) at time t = /"and Vt is the variance of the first passage time from 0 to Itt, which can be approximated by Vt — V^^J^ / {IttY. (d) Arnold 
tongue (coloured region), which shows the parameter region for synchronization to an input signal, with respect to the signal angular frequency co (vertical axis) and the signal 
strength x (horizontal axis). The dashed line is a linear approximation (equation (2.17)) of the border of the Arnold tongue when the input strength x is sufficiently small. 

integration by parts), yielding with 

2 1 f^'^ ^ 

l-^n{<p;t)=D^ni<p;t), (2.9) D = -J^ de g L/,(ef Q,(ef . (2.10) 
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Figure 4. {a) Solar radiation input of equation (2.24), where the shaded region denotes night, (b) Comparison between solar radiation input (equation (2.24)) and 
actual observed irradiance data taken from Vick & Moss [39]. The solar radiation input (solid line) and the observed data (circle) refer to left and right axes, 
respectively, (c) Gene regulatory circuit of hypothetical circadian clock. In this example, and Xj describe mRNA and protein, respectively, and x^ represses 
the transcription of Light stimulus increases the translational efficiency, (cf) Time course of Xic/ cj)) (equation (3.2)), which is a variable to be multiplied 
by the parameter p (equation (3.1)). 



See Kuramoto [18] for further details of stochastic phase 
reduction and the phase-averaging procedure. From equation 
(2.9), because n{(p = 4>-flt; t)[ = P{4>; t)] obeys a simple one- 
dimensional diffusion equation, its solution is represented by 



(2.11) 



Equation (2.11) shows that the variance of the phase after one 
period T is 

= 2DT. 

In equation (2.4), the average period corresponds to the 
mean first passage time with (/> starting from 0 to 277, and 
the period variance is the variance of the first passage time. 
Because direct calculation of the period variance is difficult, 
we approximate the period variance Vt with the phase var- 
iance V(^, after Kori et ah [26]. As the phase (/> crosses a 
threshold (/> = 277 with gradient 277/T without noise, there 
is a scaling relation y/V^ 2^y/V^T /{Itt) for sufficiently weak 
noise [26] (figure 3h). Consequently, the variance of the 
period is approximated by 



277 



d^£u-(^)2QK^)l (2.12) 



2.3. Definition of entrainability 

The entrainment property is an important characteristic 
of limit-cycle oscillators and attracts attention in systems 
biology [27-32]. For instance, a period mismatch in coupled 
oscillators is known to enhance entrainability in genetic oscil- 
lators [31]. Light stimuli affect the rate constants, i.e. the 



parameter p in equation (2.2) is perturbed as p + dp by 
the input signal. Phase d3mamics of equation (2.2) can be 
viewed as representing that of a tilted periodic potential (i.e. 
ratchet) subjected to noise. Because a synchronizable condition 
corresponds to the existence of stable points in the ratchet-like 
potential, the entrainability can be discussed without consider- 
ing the noise. Consequently, in contrast to the calculation of 
regularity, in the evaluation of the entrainability, we consider 
a case without molecular noise (i.e. (x) = 0 in equation (2.2)). 

Let be an input signal with angular frequency Con- 
sidering a weak periodic input signal dp = xpi^^)/ where x is 
the signal strength (x > 0), and applying the phase reduction 
approach to equation (2.2), the time evolution of the phase 
variable 4> is given by 



d^ 

'dt' 



N 



dxi 



,•=1 dp 
n + xZ{4>)pio)t), 



dp 



(2.13) 



with F;(<^; p) = f ,(xlc(</>); p) and 



i=l 



dp 



(2.14) 



where Z((/)) is the PRC with respect to the parameter p and 
corresponds to experimentally observed PRCs. In order to 
distinguish Z{<p) from iPRC Ui{<p), we will refer to Z{(p) as 
the parametric PRC (pPRC) [33]. Note that the definition of 
measured PRCs is different from pPRCs Z(0) in a rigorous 
definition; the experimentally measured PRCs quantify the 
phase shift Ac/) caused by light stimuli, whereas pPRCs 
Z((/)) are normalized by the strength of perturbation, i.e. 
Z((/)) = dip /dp :^A(/)/ Ap. Therefore, the ranges of the measured 



pPRCs have limitation — 7r< A(/)< tt, whereas pPRCs Z((/)) 
do not. The phase reduction can yield reliable results only 
when the perturbed trajectory is close to the unperturbed 
limit-cycle trajectory (i.e. x is sufficiently small). 

We next evaluate the extent of synchronization to the 
periodic input signal. By introducing another slow varia- 
ble \p=4>-o)t in equation (2.13), we can again apply the 
phase-averaging procedure, which yields 

dip 



:A/2 + ;^0(iA), 



(2.15) 



with Ail = fl- o) and 



277 



dez{ip+e)p{e). 



(2.16) 



The oscillator of interest can synchronize to input signals 
when there is a stable solution of ip \n if/ = 0 (equation 
(2.15)). The stable solution is an intersection point of 
and -Ail with d6)/di//<0 (an empty circle in figure 3c). 
Then, a condition for the existence of a stable solution is 



xSiirJ + n<(0< x&{>i>M) + A 



(2.17) 



where = argmax^6)(i/^) and = argmin^6)(i/^). 

We define entrainability, the extent of synchronization to 
the periodic input signal, by the width of the Arnold tongue, 
which is a domain with respect to x (signal strength) and o) 
(signal angular frequency). The shaded region in figure 3d 
represents the Arnold tongue; with parameters x ^^i^ ^ 
inside the Arnold tongue, the oscillator can synchronize to 
a periodic input signal. Because equation (2.17) constitutes 
a linear approximation of the Arnold tongue for sufficien- 
tly small Xf the width of the Arnold tongue is given by 
x{0{\pM} - ^(*Am)) under the linear approximation. Thus, we 
define the entrainability 8, or the extent of synchronization, as 

5=0(^m)-6»(«AJ. (2.18) 

Intuitively, a circadian oscillator with better entrainability 
(i.e. larger 8) can synchronize to an input signal that has a 
period further from that of the oscillator. The calculation 
above is standard in the phase reduction approach, and 
further details are available in Kuramoto [18]. 

2.4. Variational method 

We use the variational method to calculate the optimal PRCs 
which maximize the entrainability 8 subject to constant var- 
iance Vt = cr^ (the optimal solutions correspond to the 
edge in figure 1, which is described by the thick-dashed 
line). The constrained optimization of llf((/)) can be intuitively 
interpreted as maximization of weighted area (equation 
(2.18)), where the input being the weight, with constant 
area under the squared magnitude (equation (2.12)). In 
simple terms, the optimality is reached when the magnitude 
of the PRC is small during intervals when the input magni- 
tude is small (and vice versa). In the context of neuronal 
oscillators, a study [34] has used the variational method 
to calculate the optimal PRCs for stochastic synchrony 
(noise-induced synchronization [35,36]). 

The variational equation to be optimized is 

C[U]=8[U]-\Vt[U], (2.19) 

where A is the Lagrange multiplier. Note that variational 
equation (2.19) is similar to Harada et ah [37], which opti- 
mizes the input signal for the maximal entrainment under 



constant power of the input. The variational condition 
8C[U\ = 0 yields the optimal iPRC 

and the pPRC is calculated with equation (2.14): 
Z(<^) = ^E^^^^^%#^f^V- (2-21) 



Qi{<t>r 



dp 



Because i/^m and i//^ themselves depend on Ui{4>), they have 
to satisfy a self-consistent condition, i.e. equation (2.18) is 
maximal with \pM and ipm- Consequently, we maximize the 
following function: 



8{^, 8) 



77(7^ 



^(A,5), 



(2.22) 



with 



^(A, 8) 



Jo 



{p{o-A)-p{o)Y m{e+8-p) 



Qr{e + 8f 



dp 



(2.23) 



where A = ijjM- and 8 = ip^. The optimal iPRC can be 
obtained by first finding the maximum solution of ^(A, 8) 
with respect to A and 8, and then substituting the obtained 
solution \p^ = 5 and iIjm= 8-\- A into equations (2.20) and (2.21). 



2.5. Input signal of solar radiation model 

optimal PRCs depend on input signals, as seen in equations 
(2.20) and (2.21). The most common synchronizer in circadian 
oscillators is sunlight, for which the strength is determined by 
24 h-periodic solar irradiance. The solar irradiance is calcu- 
lated by 1 = Iq cos # and 1=0 when the sun is above the 
horizon (0 < tt) and below the horizon (7r< 27r), 
respectively, where # is the zenith angle and Iq is the 
maximum irradiance [38]. It can be approximated by 

p{a)t) = ramp(sin(wf)), (2.24) 

where ramp(x) is the ramp function defined by ramp(x) = x for 
X > 0 and ramp(x) = 0 for x < 0. We call equation (2.24) the solar 
radiation input, whose plot is shown in figure 4a (the shaded 
region represents night). In order to show the validity of the 
solar radiation modelling, we compare equation (2.24) with 
observed irradiance data from Vick & Moss [39], which are 
shown in a dual axis plot of figure In figure equation 
(2.24) is plotted by the soHd line (left axis) and the observed 
data by the dashed line (right axis). The solar radiation input 
of equation (2.24) is shifted horizontally, so that equation 
(2.24) becomes a good fit to the data. From figure the solar 
radiation input is in good agreement with the observed data, 
which verifies the validity of equation (2.24) as a solar radiation 
model. 

For comparison, we also use a sinusoidal input, which is 
common in nonlinear sciences: 



p{o)t) = sm{a)t). 



(2.25) 



Note that p{o)t) = B + sm{o)t), where B is an arbitrary con- 
stant, also yields the same optimal PRCs as equation (2.25), 
because a constant B in the signal is offset in equations 
(2.20) -(2.23). Although a constant B does not play any roles 
in formation of the optimal PRCs, different B result in 



different Arnold tongues in general. For calculating the 
optimal PRCs, we use equations (2.24) and (2.25). 



3. Results 

3.1. Optimal phase -response curve of solar 
radiation input 

Light stimuli generally affect the oscillatory dynamics multi- 
plicatively, i.e. they act on the rate constants or transcriptional 
efficiency of the gene regulatory circuits [3,40]. We assume 
that the ;th molecular species includes a parameter p as 

Fj{x;p) =Fj{x)+pXk, (3.1) 

where Fj{x) represents the terms that do not include p, and x^ 
is the concentration of the kih molecular species. Here, 
^E{1,2, . . . ,N} can take any value, regardless of; (both; ^ k 
and; = k are allowed). For example, let figure 4c be a gene regu- 
latory circuit of a hypothetical circadian clock, where symbols 

and H represent positive and negative regulations and Xi are 
molecular species (see Novak & Tyson [41] for typical motifs of 
biochemical oscillators). Suppose Xi and X2 are mRNA and cor- 
responding protein, respectively, and light stimuli increase the 
translational efficiency. In this case, the dynamics of light 
entrainment can be described by equation (3.1) with ; = 2, 
k=l and p being the translation rate. In equation (3.1), 
although we can also consider an alternative case 
Fj{x; p) = Fj{x) — pXk (a negative sign), the optimal pPRCs 
remain unchanged under the inversion which is seen from 
equations (2.21) and (2.23). Consequently, we consider only 
the positive case to calculate the optimal PRCs (i.e. equation 
(3.1)). However, note that relations between iPRCs and pPRCs 
are affected by the inversion of the sign, and the difference mat- 
ters when considering biological feasibility. 

When using phase reduction, the dynamics of the limit 
cycle are considered on the unperturbed limit-cycle trajec- 
tories Xlc/ and hence the points on the limit cycle can be 
uniquely determined by the phase (/>. Consequently, under 
the phase reduction, Xk is replaced by Xlc,Jc(</>) iri equation 
(3.1), where Xlc,Jc(</>) is the kih coordinate of Xlc 
(i.e. dpFj{(f)]p) = XLc,^t(</>) in equation (2.20)). Here, XLc,it(</>) 
corresponds to the time course of the concentration of the 
/cth molecular species. Because XLc,jt(</>) constitutes a core 
clock component and is generally a smooth 277-periodic 
function, we approximate it with a sinusoidal function: 

^LC,fc(</>) = 1 - 0Lsm{(j) + u), (3.2) 

where u is the initial phase and a denotes the amplitude 
of the oscillation (figure M). To ensure XLc,fc(</>) > 0/ we set 
0 < « < 1, and = 0 recovers the additive case. Because the 
initial phase u does not play any role {u is offset by 5 in equation 
(2.23)) when the white Gaussian noise is additive (i.e. Qi{x) oc 1), 
we also set w = 0. The parametric approximation of equation (3.2) 
enables an almost closed form for the overall calcu- 
lations. Although we assumed in equation (3.1) that effects of 
p only depend on Xj^, we can generalize equation (3.1) to 
Fj{x] p) = Fj{x) + pK{x), where K{x) is a nonlinear function 
and is assumed to be well approximated by 1 - asm{4> + u). 
By this generalization, our theory can be applied to other possible 
light entrainment mechanisms such as the intercellular coupling 
[42]. Our model needs only details about molecular species 
which have light input entry points but not about a whole 



molecular network. However, this advantage, in turn, means 
that we cannot specify iPRCs Ui{4>) of molecular species not 
having light input entry points. Consequently, for a noise term 
Qi{x), we assume that the white Gaussian noise is additive and 
is present only in the ;th coordinate equation (Q;((/>) = ^/q, 
where q is the noise intensity and Q/((/)) = 0 for z 7^ ;). 

Figure 5a- c shows the landscape of ^(A,5) as functions of A 
and 5, and figure 5ii-/ expresses the optimal iPRCs Uj{4>) and 
pPRCs Z( 0) for the solar radiation input (an explicit expression 
of ^(A, 8) is given in appendix A). The optimal PRC shape does 
not depend on the model parameters such as the period T, its 
variance o^, or noise intensity q. These three parameters only 
act on the magnitude of the PRCs (i.e. the vertical scaling of 
the PRCs). Consequently, we normalized T = 1, dj = \, and 
q = 1, as shown in figure 5. As the optimal PRCs depend on a, 
^(A, 8) is plotted for three cases: a = 0, (figure 5a), a = 0.5 
(figure 5h) and a: = 1.0 (figure 5c), where the maximal points 
(A, 8) yield the optimal PRCs using equations (2.20) and 
(2.21). The maximal parameters A and 8 are calculated numeri- 
cally. Figure 5d-f describes the optimal iPRCs (solid line) and 
pPRCs (dashed line) for a = 0, 0.5 and 1.0, respectively. When 
a = 0, i.e. the input signal is additive, ^(A,5) achieves a maxi- 
mum for A = 77 and arbitrary 8, yielding sinusoidal PRCs as 
the optimal solution (figure 5d). Although the input signal 
p{4>)^s riot sinusoidal, the optimal PRCs obtained using the vari- 
ational method become sinusoidal. In other words, considering 
optimality, resonator-type oscillators have an advantage over 
integrator-t3rpe oscillators. For a>0, the input signal p{4>) 
depends on the concentration of the /cth molecular species. 
From figure 5h, the optimal parameters for o: = 0.5 are (A, 8) = 
(2.31, 1.99) and (3.98, 4.30), which are different from A = 77 
(these two sets yield symmetric PRCs with respect to the hori- 
zontal axis). Figure 5e shows the optimal iPRCs Uj{4>) and 
pPRCs Z((/)) for q: = 0.5. Interestingly, the optimal iPRCs 
and pPRCs for = 0.5 have a dead zone (region of 1 ^ (/> ^ 2 
in figure 5e) in which the input signal neither advances nor 
delays the clock. From equations (2.20) - (2.21 ) and the solar radi- 
ation input of equation (2.24), the optimal PRCs inevitably 
include a dead zone if the optimal A is not tt. For a= 1.0, 
there are four sets of parameters (A, 8) that give optimal PRCs: 
(2.30, 2.72), (2.30, 1.26), (3.98, 3.56) and (3.98, 5.02) (PRCs with 
these four sets are symmetric with each other with respect to 
the horizontal axis or (/) = 37r/2). Consequently, the optimal 
PRCs shown in figure 5/ have a dead zone as in the case of 
a 0.5. 

3.2. Dead zone length 

From the results discussed above, the optimal PRCs have a 
dead zone when > 0. We next studied the length of the 
dead zone as a function of a (figure 6a) and improvements 
in the entrainability induced by the dead zone (figure 6h) 
for the solar radiation input. Because the dead zone, which 
is a null interval in PRCs, emerges when the optimal 
parameter is A 7^ 77, we can naturally define its length as 

L = \\-M, (3.3) 

where A is the maximum value of ^(A, 8). As seen in 
figure 6a, a dead zone clearly exists when a> 0, and the 
length increases with increasing a for < 0.8. Even for 
= 0.1, when the oscillation amplitude of XLc,fc(</>) (the con- 
centration of a molecular species modulated by the light- 
sensitive parameter p; cf. figure M) is very small, we observe 




4^ 
b 



OO 



Figure 5. (a-c) Landscape of ^(A, 8) as a function of A and 8 with solar radiation input (equation (2.24)) for (a) a = 0, (b) a = 0.5 and (c) a = 1.0, where 
the maximum points are parameters for the optimal PRC. (d-f) Optimal PRCs with solar radiation input: (d) a = 0, (e) a = 0.5 and (f) « = 1.0. In (d-f), 
the solid and dashed lines denote iPRCs Ujicf)) and pPRCsZ((/)), respectively (in (d), solid and dashed lines are indistinguishable). The maximal parameters (A,8) 
for (d-f) are (d) (tt, 0), (e) (2.31, 1.99) and (f) (2.30, 2.72). In [d], a parallel shift of the PRC is also optimal (8 can be an arbitrary value). In (e), symmetric 
PRCs with respect to the horizontal axis are also optimal. In (f), symmetric PRCs with respect to the horizontal axis or (/> = 37r/2 are also optimal (see the text). 
The pPRCs correspond to experimentally observed PRCs. 




a dead zone with a length of L = 0.475, which corresponds 
to about 3 h within 24 h, indicating the universality of 
having a dead zone in order to attain optimality. The 
improvement in the entrainability that is induced by a dead 
zone is calculated by comparing the entrainability of the 
optimal PRCs with that of typical sinusoidal PRCs. We con- 
sider sinusoidal functions for both the iPRC Uj{(f)) and 
pPRC Z(0) by setting 

Uy(c/>) ocsin(c/> + c) (3.4) 

and 

Z((/>) ocsin((/> + c), (3.5) 

where c is the parameter to be optimized so that entrainabil- 
ity is maximized for each a (see appendix B for the explicit 
expressions). Equations (3.4) and (3.5) are scaled, so that 



they satisfy the constraints on the period variance (equation 
(2.12)). We calculated the ratios 

Ru = ^ and Rz=^. (3.6) 

where 8u and 8z represent the entrainabilities for the cases of 
the sinusoidal iPRC and pPRC, respectively, calculated 
for the solar radiation input. For the sinusoidal iPRC of 
equation (3.4), the entrainability is calculated with pPRC 
via equation (2.14). Ru and Rz quantify the improvement 
rate of the optimal PRCs over the sinusoidal iPRC {Ru) and 
pPRC {Rz)- In figure 6h, the solid and dashed lines show 
Ru and Rz, respectively, as a function of a. Both ratios mono- 
tonically increase as a increases, which shows that the 
optimal PRC with a dead zone exhibits better entrainability 
when the oscillation of Xlc,Jc(</>) has a larger amplitude. 
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Figure 7. (a-c) Landscape of ^(A, 8) as functions of A and 8 with sinusoidal input for (a) a = 0, (b) a = 0.5 and (c) a = 1.0, where the maximum points 
are parameters for the optimal PRC. (d-f) Optimal PRCs with sinusoidal input (equation (2.25)): (d) a = 0, (e) a = 0.5 and (f) « = 1.0, where the solid and 
dashed lines denote iPRCs Ujicf)) and pPRCsZ((/)). The maximal parameter (A, 8) is (tt, 0) in all cases. In (d), a parallel shift of the PRC is also optimal (8 can be 
an arbitrary value). PRCs that are symmetric with respect to the horizontal axis are also optimal. The pPRCs correspond to experimentally observed PRCs. 



When the concentration of Xi^cA4>) is low, the effects of the 
input signal on the circadian oscillators are smaller. This is 
because pPRC Z{4>), which quantifies the extent of the 
phase shift owing to the stimulation of the parameter, 
depends on the concentration XLc,fc(</>) (see equation (2.14)). 
However, even within the range cf) where Xi^c,k{^) ^^s smaller 
values, the iPRC Uj{4>) contributes to an increase in the var- 
iance of the period, regardless of the concentration. From this, 
we see that having an iPRC with a smaller magnitude when 
the concentration of Xi^c,k{^) is smaller results in a smaller 
variance, which results in a larger entrainability for a constant 
variance of the period. Although this qualitatively explains 
the benefit of a dead zone, for some input values, the optimal 
PRCs may not contain a dead zone for any value of a. This 
will be shown in the following. 

3.3. Optimal phase -response curve of sinusoidal input 

Because the optimal PRCs depend on input signals (equations 

(2.20) and (2.21)), we next consider a typical periodic input 
signal, a sinusoidal function (equation (2.25)). In this case, 
^(A, 8) is calculated in a closed form (an explicit expression 
of ^(A, 8) is given in appendix A), which is plotted as functions 
of A and 8 in figure 7a-c for three cases: a = 0 (figure 7a), a = 
0.5 (figure 7h) and = 1.0 (figure 7c). As can been seen from 
figure 7a- c, ^(A, 8) yields the maximal value for (A, 8) = {tt, 
utt) for 0 < < 1, where n is an integer and when a = 0, 8 
can take any value. Figure 7d-f expresses the optimal iPRCs 
Uj{4>) and pPRCs Z{4>) for the sinusoidal input. For a = 0, 
the optimal PRC is sinusoidal (figure 7d) and for = 0.5, the 
optimal PRC is still close to a sinusoidal function (figure 7e). 
When increasing ato a= 1.0, the PRC diverges from the sinu- 
soidal function and exhibits almost positive values (figure 7/ ). 
We see that the optimal PRCs owing to equations (2.20) and 

(2.21) do not exhibit a dead zone for any o^-values (figure 
7d-f) when the input signal is a simple sinusoidal function. 



4. Discussion 

The existence of a dead zone optimizes both entrainability and 
regularity. It is rather obvious that optimization of regularity 
alone leads to a dead zone [43], because null response 
means no effect by any kind of fluctuations. Our result instead 
shows that optimality of both entrainability and regularity, 
which are in a trade-off relationship, is uniquely achieved by 
a dead zone. Our finding is fairly general, because a dead 
zone always exists in an optimal PRC unless a = 0 (additive 
stimulation). Along with the fact that T, ctt and q affect only 
the scaling of the optimal PRCs, when the input signal affects 
the dynamics multiplicatively (i.e. a > 0), the existence of a 
dead zone always provides a synchronization advantage. 
This is supported by many experimental studies of various 
species that report the existence of a dead zone in the PRC 
[1] (cf. figure 2c, d). Our general result suggests that circadian 
oscillators have fully adapted to solar radiation to improve 
synchronization. Indeed, many experimental findings imply 
that circadian oscillators have adapted to actual solar radiation 
[44]: for various animals, light -dark (LD) cycles that include a 
twilight period result in better entrainability than do abrupt 
LD cycles (on-off protocols) [44]. In this regard, another inter- 
esting problem is optimal entrainment [37] of circadian clocks 
by light stimuli. As two different input signals, the solar radi- 
ation and sinusoidal inputs, yield the same optimal PRCs for 
a = 0, optimal inputs and optimal PRCs do not have one-to- 
one correspondence. Thus, the optimal inputs are not trivial 
and this problem should be pursued in our future studies. 

The solar radiation input plays an essential role, because it 
yields a dead zone in the optimal PRC, whereas a sinusoidal 
signal does not (figure 7). In other words, oscillators that are 
entrained by stimuli other than solar radiation may not exhibit 
a dead zone in their PRCs. This is indeed found in mammals. 
Mammals possess a master clock in their suprachiasmatic 
nucleus (SCN), which receives light stimuli via retinal 



photoreceptors, and peripheral clocks in body cells [45]. The 
peripheral oscillators are entrained by several stimuli such as 
feeding and signals from the SCN through chemical pathways 
(e.g. hormones) [45,46]. By injection experiments of a hormone, 
Balsalobre et al. [47] reported that the PRCs of the peripheral 
oscillators in the liver do not have a dead zone. 

Our result also agrees with other experimental observations. 
Our theory implies that a dead zone should be located where the 
concentration XLc,fc( 4>) is low (0 < (/> < tt in figure M), and that to 
achieve optimality, the concentration of Xi^c,k{4>) should be maxi- 
mal in the region where the PRCs exhibit a large phase shift. In 
Drosophila, the timeless (tim) gene is regarded as the molecular 
implementation of Xi^c,k{4>)' It is experimentally known that 
light enhances the degradation of the gene product (the TIM 
protein) [48,49], and the TIM protein peaks during the late eve- 
ning. Figure 2c shows observations of the PRC of Drosophila 
against light pulses as a function time from Hall & Rosbash 
[19]; circles describe the experimental data, and the solid line 
expresses a trigonometric curve fitting (fourth order). Because 
the centre of the part of the PRC that can be phase shifted approxi- 
mately corresponds to the peak of the concentration, as denoted 
above, when estimated from the PRC alone, the concentration 
peak of the TIM protein should occur at about 18 h. This time is 
also close to the experimental evidence (i.e. late evening). There- 
fore, our theory can be used to h3^othesize further molecular 
behaviour affected by light stimuli. 

In summary, we have constructed a model that regards 
circadian oscillators as a global optimization of entrainability 
and regularity. We have shown that our model is consistent 
with much experimental evidence as mentioned above. The 
extension and improvement of our method are possible and 
they are left as an area of future study. 
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Appendix A, Explicit expression of ^A, 8) 
A.I. Solar radiation input case 

For the solar radiation input case (equation (2.24)), ^(A, 8) is 
given by 



(Al) 



with 



%{^, 8) =-r^[-3a^ sin(3A + 25) + 32c^cos(2A + 8) 

+12a^A cos(25 + A) - 6a^ sin(2§ + A) 

+ 120:^77 cos (A + 8f - 128q:cos(A + 8) 

+ I— 24q^ TTCOS (5) V (128c^ + 18c^2 sin 5) cos o 

+(-1277 + 24A)c^2 - 4877+ 48A} cos A 

+ 12 Q sin A + 77^ cos {8f 

+ (24c^^ 77sin A sin 5 - 32a) cos 8 

+f— 64a sin 5 - 48 - 27q^^) sin a + 4877 + 12^^77], 



where ^^(A, 8) = + 277, -8 + 277). We show equation 

(A 1) as functions of A and 8 in figure 5d-f. 



A.2, Sinusoidal input case 



For the sinusoidal input case (equation (2.25)), ^(A, 8) is 
given by 

^(A, 8)=^(1- cos ^)[-a^ cos(25 + A) + 2a^ + 4]. (A 2) 

We plot equation (A 2) as functions of A and 8 in 
figure 7d-f. 



Appendix B. Explicit expression of sinusoidal 
phase -response curves 

B.I. Sinusoidal infinitesimal phase -response 
curve 

An explicit expression for the sinusoidal iPRC (equation (3.4)) is 



which yields the period variance of Vj = d^. Then, the 
corresponding pPRC is given by 



Z(</.) = J^^sin(<^ + c)(l-asin</.), (B2) 



where we used equation (2.14). 



B.2. Sinusoidal parametric phase -response 
curve 

For the pPRC Z{4>) to be a sinusoidal function, the iPRC 
Uj{(p) must be 



dp 



sin{4> + c) 



(B3) 



where we used equation (2.14). An explicit expression of 
equation (B 3) is 

1 sin((/) + c) 



where J\f{c) is a normalizing term 



(B4) 



_ 1 Vl -2Vl - c^^ +2| cos(2c) 



,2^3/2 



Equation (B4) is normalized, so that the period variance 
becomes Vj = d^. Using equation (2.14), the corresponding 
pPRC is a sinusoidal function: 



2(c/>) =-^sin((/> + c) 



(B5) 



which is an explicit expression of the sinusoidal pPRC 
(equation (3.5)). 
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